tidymodels
El objetivo de un modelo es proporcionar un resumen simple de baja dimensión de un conjunto de datos. Idealmente, el modelo capturará patrones generados por el fenómeno de interés e ignorará el “ruido” (es decir, la variación aleatoria que no le interesa).
Tradicionalmente, el enfoque del modelado está en la inferencia o confirmar que una hipótesis es cierta. Hay un par de ideas que nos pueden ayudar a hacer inferencias correctamente:
Cada observación puede usarse para exploración o confirmación, no ambas.
Puede utilizar una observación tantas veces como desee para la exploración, pero solo puede utilizarla una vez para la confirmación. Tan pronto como utilice una observación dos veces, pasará de la confirmación a la exploración.
Esto es necesario porque para confirmar una hipótesis debe utilizar datos independientes de los datos que utilizó para generar la hipótesis. De lo contrario, será demasiado optimista. No hay absolutamente nada de malo en la exploración, pero nunca debe vender un análisis exploratorio como un análisis confirmatorio porque es fundamentalmente engañoso.
Si realmente desea realizar un análisis confirmatorio, un enfoque es dividir sus datos en tres partes antes de comenzar el análisis:
El 60% de sus datos se destina a un conjunto de entrenamiento, training, (o exploración). Puede hacer lo que quiera con estos datos: visualizarlos y ajustarles toneladas de modelos.
El 20% se destina a un conjunto de consultas, query. Puede usar estos datos para comparar modelos o visualizaciones a mano, pero no puede usarlos como parte de un proceso automatizado.
El 20% se retiene para una prueba testing. Solo puede usar estos datos UNA VEZ, para probar su modelo final.
Esta partición le permite explorar los datos de entrenamiento, generando ocasionalmente hipótesis candidatas que verifica con el conjunto de consultas. Cuando esté seguro de que tiene el modelo correcto, puede verificarlo una vez con los datos de prueba.
Tidymodels, es una interfaz que unifica bajo un único marco cientos de funciones de distintos paquetes, facilitando en gran medida todas las etapas de preprocesado, entrenamiento, optimización y validación de modelos predictivos. Los paquetes principales que forman parte del ecosistema tidymodels son:
parsnip: para la definición de modelos.
Implementación tidy del caret.
recipes: para el preprocesado de datos y feature engineering.
rsample: para validar los modelos por métodos de resampling.
dials: para crear y manejar el valor de los hiperparámetros.
tune: para hacer tuning de modelos.
yardstick: para calcular métricas de modelos.
workflows: para combinar todos los pasos del preprocesado y modelado en un único objeto.
El procesamiento lo realizan los paquetes rsample y
recipes. El modelo principal está en parsnip
(el equivalente de Tidy a caret) y para la
validación Yardstick. Cabe señalar que caret
sigue teniendo la mejor implementación de las tablas de matrices de
confusión, con un paquete adaptado
ConfusionTableR.
Evaluar la capacidad predictiva de un modelo consiste en comprobar cómo se acercan las predicciones a los valores verdaderos de la variable respuesta. Para poder cuantificarlo de forma correcta, se necesita disponer de un conjunto de observaciones, de las que se conozca la variable respuesta, pero que el modelo no haya “visto”, es decir, que no hayan participado en su ajuste. Con esta finalidad, se divide el conjunto de datos disponibles en un conjunto de entrenamiento (train) y un conjunto de prueba (test).
El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 80%-20% suele dar buenos resultados. El reparto debe hacerse de forma aleatoria o aleatoria-estratificada. Evaluación de modelos (ver link).
Para este ejemplo, analizaremos el conjunto de datos
SaratogaHiyses, que se encuentra en el paquete
nisaucData, el cual contiene información sobre los precio
de 1728 viviendas de Saratoga County, New York, USA en el año 2006.
Además del precio, este dataset incluye otras 15 variables:
price: precio de la vivienda.
lotSize: metros cuadrados de la vivienda.
age: antigüedad de la vivienda.
landValue: valor del terreno.
livingArea: metros cuadrados habitables.
pctCollege: porcentaje del vecindario con título universitario.
bedrooms: número de dormitorios.
firplaces: número de chimeneas.
bathrooms: número de cuartos de baño (el valor 0.5 hace referencia a cuartos de baño sin ducha).
rooms: número de habitaciones.
heating: tipo de calefacción.
fuel: tipo de alimentación de la calefacción (gas, electricidad o diesel).
sewer: tipo de desagüe.
waterfront: si la vivienda tiene vistas al lago.
newConstruction: si la vivienda es de nueva construcción.
centralAir: si la vivienda tiene aire acondicionado.
El objetivo es obtener un modelo capaz de predecir el precio del alquiler casa. Por lo tanto para poder llamar los datos tomamos el siguiente código:
data("SaratogaHouses", package = "mosaicData")
data <- SaratogaHouses
colnames(data) [1] "price" "lotSize" "age" "landValue" "livingArea"
[6] "pctCollege" "bedrooms" "fireplaces" "bathrooms" "rooms"
[11] "heating" "fuel" "sewer" "waterfront" "newConstruction"
[16] "centralAir"
Con este conjunto de datos procedemos a hacer la división entre train y test.
library(tidymodels)Registered S3 method overwritten by 'tune':
method from
required_pkgs.model_spec parsnip
── Attaching packages ────────────────────────────────────────────────────────────── tidymodels 0.1.3 ──
✓ broom 0.7.8 ✓ recipes 0.1.16
✓ dials 0.0.9 ✓ rsample 0.1.0
✓ dplyr 1.0.7 ✓ tibble 3.1.2
✓ ggplot2 3.3.5 ✓ tidyr 1.1.3
✓ infer 0.5.4 ✓ tune 0.1.5
✓ modeldata 0.1.0 ✓ workflows 0.2.2
✓ parsnip 0.1.6 ✓ workflowsets 0.0.2
✓ purrr 0.3.4 ✓ yardstick 0.0.8
── Conflicts ───────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
x purrr::discard() masks scales::discard()
x dplyr::filter() masks stats::filter()
x dplyr::lag() masks stats::lag()
x recipes::step() masks stats::step()
• Use tidymodels_prefer() to resolve common conflicts.
# Reparto de datos en train y test
set.seed(123)
split_inicial <- initial_split(
data = data, # datos que queremos dividir
prop = 0.8, # Proporción del 80%
strata = price # Variable respuesta
)
data_train <- training(split_inicial)
data_test <- testing(split_inicial)Es importante verificar que la distribución de la variable respuesta
es similar en el conjunto de entrenamiento y en el de test. Para
asegurar que esto se cumple, la función initial_split()
permite identificar con el argumento strata la variable en
base a la cual hacer el reparto.
summary(data_train$price) Min. 1st Qu. Median Mean 3rd Qu. Max.
5000 145000 189050 212423 259000 775000
summary(data_test$price) Min. 1st Qu. Median Mean 3rd Qu. Max.
20000 145000 189950 210157 257750 620000
Este tipo de reparto estratificado asegura que el conjunto de entrenamiento y el de test sean similares en cuanto a la variable respuesta, sin embargo, no garantiza que ocurra lo mismo con los predictores. Por ejemplo, en un set de datos con 100 observaciones, un predictor binario que tenga 90 observaciones de un grupo y solo 10 de otro, tiene un alto riesgo de que, en alguna de las particiones, el grupo minoritario no tenga representantes. Si esto ocurre en el conjunto de entrenamiento, algunos algoritmos darán error al aplicarlos al conjunto de test, ya que no entenderán el valor que se les está pasando. Este problema puede evitarse eliminando variables con varianza próxima a cero, lo cual veremos más adelante.
No se deben incluir en el modelo predictores que contengan un único
valor (varianza-cero) ya que no aportan información. Tampoco es
conveniente incluir predictores que tengan una varianza próxima a cero,
es decir, predictores que toman solo unos pocos valores, de los cuales,
algunos aparecen con muy poca frecuencia. El problema con estos últimos
es que pueden convertirse en predictores con varianza cero cuando se
dividen las observaciones por validación cruzada o
bootstrap.
¿Recuerdas que lo vimos en el análisis exploratorio, EDA?
Aquí, también podemos utilizar la función step_nzv() del
paquete recipes que identifica predictores potencialmente
problemáticos, es decir, aquellos que tiene un único valor o varianza-
cero que cumplen dos condiciones:
Variación de frecuencia: Es la variación entre la frecuencia del
valor más común y la frecuencia del según valor más común. Este ratio
tiende a 1 si las frecuencias están equidistribuidas y a valores grandes
cuando la refuencia del valor mayoritario supera por mucho al resto (el
denominador es un número decimal pequeño). Valor por defecto
freq_cut = 95/5.
Porcentaje de valores únicos: número de valores únicos divido
entre el total de la muestra (multiplicado por 100). Esta valor se
acerca a cero cuando la varianza-cero es mayor. Valor por defecto
uniqueCut = 10.
Si bien, cuando eliminamos predictores no informativos podría considerarse un paso propio del proceso de selección de predictores, dado que consiste en un filtrado por varianza, tiene que realizarse antes de estandarizar los datos, ya que después, todos los predictores tienen varianza 1.
Cuando los predictores son numéricos, la escala en la que se miden, así como la magnitud de su varianza pueden influir en gran medida en el modelo. Muchos algoritmos de machine learning (SVM, redes neuronales, lasso, etc) son sensibles a esto, de forma que, si no se igualan de alguna forma los predictores, aquellos que se midan en una escala mayor o que tengan más varianza dominarán el modelo aunque no sean los que más relación tienen con la variable respuesta. Existen principalmente 2 estrategias para evitarlo:
Centrado: consiste en restarle a cada valor la media del predictor al que pertenece. Si los datos están almacenados en un dataframe, el centrado se consigue restándole a cada valor la media de la columna en la que se encuentra. Como resultado de esta transformación, todos los predictores pasan a tener una media de cero, es decir, los valores se centran en torno al origen.
Normalización (estandarización): consiste en transformar los datos de forma que todos los predictores estén aproximadamente en la misma escala. Hay dos formas de lograrlo:
\[ z = \frac{x - \mu}{\sigma} \]
\[ x_{nom}=\frac{x-x_{min}}{x_{max}-x_{min}} \] NOTA: Nunca se debe estandarizar las variables después de ser binarizadas.
La binarización consiste en crear nuevas variables dummy con
cada uno de los niveles de las variables cualitativas. A este proceso
también se le conoce como one hot encoding. Por ejemplo, una
variable llamada color que contenga los niveles rojo, verde y azul, se
convertirá en tres nuevas variables (color_rojo,
color_verde, color_azul), todas con el valor 0
excepto la que coincide con la observación, que toma el valor 1.
Por defecto, la función
step_dummy(all_nominal(), -all_outcomes()) binariza todas
las variables almacenadas como tipo factor o
character, excepto la variable respuesta. Además, elimina
uno de los niveles para evitar redundancias. Volviendo al ejemplo
anterior, no es necesario almacenar las tres variables, ya que, si
color_rojo y color_verde toman el valor 0, la variable color_azul toma
necesariamente el valor 1. Si color_rojo o
color_verde toman el valor 1, entonces color_azul es
necesariamente 0.
El paquete recipes incorpora una amplia variedad de
funciones para preprocesar los datos, facilitando el aprendizaje de las
transformaciones únicamente con observaciones de entrenamiento, y poder
aplicarlas después a cualquier conjunto de datos. La idea detrás de este
paquete es la siguiente:
Definir cuál es la variable respuesta, los predictores y el
conjunto de datos de entrenamiento, recipe().
Definir todas las transformaciones (escalado, selección,
filtrado…) que se desea aplicar, step_().
Aprender los parámetros necesarios para dichas transformaciones
con las observaciones de entrenamiento rep().
Aplicar las transformaciones aprendidas a cualquier conjunto de
datos juice(), bake().
# Se almacenan en un objeto `recipe` todos los pasos de preprocesado y, finalmente, se aplican a los datos.
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
transformerData Recipe
Inputs:
Operations:
Removing rows with NA values in all_predictors()
Sparse, unbalanced variable filter on all_predictors()
Centering for all_numeric(), -all_outcomes()
Scaling for all_numeric(), -all_outcomes()
Dummy variables from all_nominal(), -all_outcomes()
Una vez que se ha definido el objeto recipe, con la función
prep() se aprenden las transformaciones con los datos de
entrenamiento y se aplican a los dos conjuntos con
bake().
# Se entrena el objeto recipe
transformer_fit <- prep(transformer)
# Se aplican las transformaciones al conjunto de entrenamiento y de test
data_train_prep <- bake(transformer_fit, new_data = data_train)
data_test_prep <- bake(transformer_fit, new_data = data_test)
glimpse(data_train_prep)Rows: 1,380
Columns: 18
$ lotSize <dbl> -0.43894231, 0.26607619, 0.48189818, 0.61139137, 12.19383812, -0.51088…
$ age <dbl> 3.62268073, 0.11471419, 0.28667334, -0.91704067, -0.50433872, -0.64190…
$ landValue <dbl> -0.7674967, -0.5829431, -0.3570715, -0.3543170, -0.8363600, -0.9603140…
$ livingArea <dbl> 0.31203819, -0.96981699, -0.19293506, -0.20588309, -1.69490679, -0.730…
$ pctCollege <dbl> -0.4558859, -3.2928560, -0.4558859, -0.4558859, -1.4341515, -1.4341515…
$ bedrooms <dbl> 1.0465989, 1.0465989, -0.1681248, -0.1681248, -1.3828485, -0.1681248, …
$ fireplaces <dbl> 0.7295147, 0.7295147, -1.0876880, -1.0876880, -1.0876880, -1.0876880, …
$ bathrooms <dbl> -1.3602353, -1.3602353, -0.6045490, 0.1511373, -1.3602353, -0.6045490,…
$ rooms <dbl> 0.425229547, 0.425229547, 0.425229547, -0.440920675, -1.307070897, 0.4…
$ price <int> 109000, 120000, 90000, 120000, 85860, 127000, 89900, 60000, 87500, 112…
$ heating_hot.water.steam <dbl> 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1…
$ heating_electric <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ fuel_electric <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ fuel_oil <dbl> 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ sewer_public.commercial <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1…
$ sewer_none <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ newConstruction_No <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
$ centralAir_No <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1…
Tras el preprocesado de los datos, se han generado un total de 18 variables (17 predictores y la variable respuesta).
Ya analizado los datos (entenderlos por medio del EDA), que han sido preprocesados y los predictores seleccionados, el siguiente paso es emplear un algoritmo de machine learning que permita crear un modelo capaz de representar los patrones presentes en los datos de entrenamiento y generalizarlos a nuevas observaciones. Encontrar el mejor modelo no es fácil, existen multitud de algoritmos, cada uno con unas características propias y con distintos parámetros que deben ser ajustados. Por lo general, las etapas seguidas para obtener un buen modelo son:
Ajuste/entrenamiento: consiste en aplicar un algoritmo de machine learning a los datos de entrenamiento para que el modelo aprenda.
Evaluación/validación: el objetivo de un modelo
predictivo no es ser capaz de predecir observaciones que ya se conocen,
sino nuevas observaciones que el modelo no ha visto. Para poder estimar
el error que comete un modelo es necesario recurrir a estrategias de
validación entre las que destacan: validación simple,
bootstrap y validación cruzada.
Optimización de hiperparámetros: muchos
algoritmos de machine learning contienen en sus ecuaciones uno o varios
parámetros que no se aprenden con los datos, a estos se les conoce como
hiperparámetros. Por ejemplo, los SVM
lineal tiene el hiperparámetro de coste C y los árboles
de regresión la profundidad del árbol tree_depth y el
número mínimo de observaciones por nodo min_n. No existe
forma de conocer de antemano cuál es el valor exacto de un
hiperparámetro que da lugar al mejor modelo, por lo que se tiene que
recurrir a estrategias de validación para comparar distintos
valores.
Predicción: una vez creado el modelo, este se emplea para predecir nuevas observaciones.
Es a lo largo de todo este proceso donde más destacan las
funcionalidades ofrecidas por tidymodels, permitiendo
emplear la misma sintaxis para ajustar, optimizar, evaluar y predecir un
amplio abanico de modelos, variando únicamente el nombre del
algoritmo.
Aunque tidymodels permite todo esto con apenas unas
pocas líneas de código, son muchos los argumentos que pueden ser
adaptados, cada uno con múltiples posibilidades. Con el objetivo de
exponer mejor cada una de las opciones, en lugar de crear directamente
un modelo final, se muestran ejemplos de menor a mayor complejidad.
Los modelos de tidymodels son accesibles por medio del
paquete parsnip. Aunque un mismo modelo puede estar en
varios paquetes con nombres distintos.
Para crear un modelo se requieren únicamente tres pasos:
Definir el tipo de modelo y sus parámetros
Definir qué implementación del modelo se quiere aplicar
(engine)
Ajustar el modelo
Vamos a ajustar un modelo de árbol
de regresión para predecir el precio de la vivienda de Saratoga en
función de todos los predictores disponibles.
A excepción de la implementación del algoritmo ("rpart"),
rpart (viene por default) o "C5.0" (solo para
clasificación), todos los demás argumentos de la función
decision_tree() se dejan por defecto.
model_tree <- decision_tree(mode = "regression") %>%
set_engine(engine = "rpart")
model_treeDecision Tree Model Specification (regression)
Computational engine: rpart
Con las líneas anteriores sólo guardamos el tipo de algoritmo (sus parámetros, hiperparámetros y el paquete en el que está implementado). El sigiente paso es ajustar el modelo, para los cual tenemos dos opciones:
fit() que emplea como argumento una formula para
definir la variable respuesta y los predictores.
fit_xy() que emplea los argumentos x y
y para definir la matriz de predictores y el vector con la
variable respuesta.
# Entrenamiento empleando fórmula
model_tree_fit <- model_tree %>%
fit(
formula = price ~ .,
data = data_train_prep
)# Entrenamiento empleando x e Y.
variable_respuesta <- "price"
predicores <- setdiff(colnames(data_train_prep), variable_respuesta)
model_tree_fit <- model_tree %>%
fit_xy(
x = data_train_prep[, predicores],
y = data_train_prep[[variable_respuesta]]
)Para ver nuestro modelo utilizamos lo que se almacena dentro de
model_tree_fit como fit:
model_tree_fit$fitn= 1380
node), split, n, deviance, yval
* denotes terminal node
1) root 1380 1.382982e+13 212423.0
2) livingArea< 0.3856801 956 3.951824e+12 172884.3
4) landValue< -0.3405443 540 1.320059e+12 147526.2
8) livingArea< -0.3968666 397 5.752973e+11 135085.3 *
9) livingArea>=-0.3968666 143 5.127270e+11 182065.0 *
5) landValue>=-0.3405443 416 1.833788e+12 205801.0
10) landValue< 0.8136044 359 1.118047e+12 196236.6
20) livingArea< -0.001951587 251 7.002854e+11 182517.3 *
21) livingArea>=-0.001951587 108 2.607222e+11 228121.3 *
11) landValue>=0.8136044 57 4.760637e+11 266039.8 *
3) livingArea>=0.3856801 424 5.013749e+12 301571.5
6) landValue< 0.9733672 311 1.944443e+12 268743.7
12) livingArea< 1.830195 277 1.304977e+12 256947.8
24) centralAir_No>=0.5 140 4.264425e+11 233996.8 *
25) centralAir_No< 0.5 137 7.294296e+11 280401.5 *
13) livingArea>=1.830195 34 2.869182e+11 364844.9 *
7) landValue>=0.9733672 113 1.811733e+12 391920.8
14) rooms< 1.940992 89 1.063606e+12 365619.8
28) age< -0.8998448 41 9.079364e+10 317584.5 *
29) age>=-0.8998448 48 7.974026e+11 406650.0
58) livingArea< 0.7029069 13 7.632773e+10 320607.7 *
59) livingArea>=0.7029069 35 5.890850e+11 438608.6
118) age>=-0.7622774 27 2.569939e+11 403566.7 *
119) age< -0.7622774 8 1.870413e+11 556875.0 *
15) rooms>=1.940992 24 4.582587e+11 489453.6 *
Lo que buscamos ver es que tan certera es la predicción de nuestro modelo, ya sea con observaciones que no ha visto el modelo y por ende, con observaciones futuras.
El error mostrado por defecto tras entrenar un modelo suele ser el error de entrenamiento, el error que comete el modelo al predecir las observaciones que ya ha “visto”.
Estos errores son útiles para entender cómo está aprendiendo el modelo (estudio de residuos), no es una estimación realista de cómo se comporta el modelo ante nuevas observaciones (el error de entrenamiento suele ser demasiado optimista).
Para conseguir una estimación más certera, y antes de recurrir al
conjunto de test, se pueden emplear estrategias de validación
basadas en resampling. El paquete rsampler
incorpora los métodos Bootstrap (bootstraps), V-Fold
Cross-Validation y Repeated V-Fold Cross-Validation
(vfold_cv), Nested or Double Resampling
(nested_cv), Group V-Fold Cross-Validation
(group_vfold_cv), Leave-One-Out Cross-Validation
(loo_cv), Monte Carlo Cross-Validation(mc_cv).
Cada uno funciona internamente de forma distinta, pero todos ellos se
basan en la idea: ajustar y evaluar el modelo de forma repetida,
empleando cada vez distintos subconjuntos creados a partir de los datos
de entrenamiento y obteniendo en cada repetición una estimación del
error (simulaciones de los datos). El promedio de todas las
estimaciones tiende a converger en el valor real del error de test.
Se ajusta de nuevo el modelo, esta vez con validación cruzada repetida para estimar su error.
En primer lugar, se crea un objeto resample que contenga
la información sobre las observaciones que forman parte de cada
partición. Dado que el reparto es aleatorio, es importante emplear una
semilla set.seed() para que los resultados sean
reproducibles.
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
repeats = 10,
strata = price
)
head(cv_folds)Con la función fit_resamples() el modelo se ajusta y
evalúa con cada una de las particiones de forma automática, calculando y
almacenando en cada iteración la métrica de interés. Como se comentó
anteriormente, cuando se realizan validaciones por resampling,
el preprocesado debe ocurrir dentro de cada iteración.
Para seguir este principio, las particiones se crean con el conjunto de datos de entrenamiento sin preprocesar y se le pasa al argumento preprocessor un objeto recipe que contenga la definición del preprocesado.
model_tree <- decision_tree(mode = "regression") %>%
set_engine(engine = "rpart")
validation_fit <- fit_resamples(
object = model_tree,
# El objeto recipe no tiene que estar entrenado
preprocessor = transformer,
# Las resamples se tienen que haber creado con los datos sin
# prerocesar
resamples = cv_folds,
metrics = metric_set(rmse, mae),
control = control_resamples(save_pred = TRUE)
)
head(validation_fit)Los resultados de fit_resamples() se almacenan en forma
de tibble,
donde las columnas contienen la información sobre cada partición: su id,
las observaciones que forman parte, las métricas calculadas, si ha
habido algún error o warning durante el ajuste, y las predicciones de
validación si se ha indicado
control = control_resamples(save_pred = TRUE). La
información puede extraerse haciendo un unnest() de las
columnas o empleando las funciones auxiliares
collect_predictions() y collect_metrics().
# Métricas promedio de todas las particiones
validation_fit %>% collect_metrics(summarize = TRUE)# Métricas individuales de cada una de las particiones
validation_fit %>% collect_metrics(summarize = FALSE) %>%
head()library("ggpubr")Registered S3 method overwritten by 'data.table':
method from
print.data.table
# Valores de validación (mae y rmse) obtenidos en cada partición y repetición.
p1 <- ggplot(
data = validation_fit %>% collect_metrics(summarize = FALSE),
aes(x = .estimate, fill = .metric)) +
geom_density(alpha = 0.5) +
theme_bw()
p2 <- ggplot(
data = validation_fit %>% collect_metrics(summarize = FALSE),
aes(x = .metric, y = .estimate, fill = .metric, color = .metric)) +
geom_boxplot(outlier.shape = NA, alpha = 0.1) +
geom_jitter(width = 0.05, alpha = 0.3) +
coord_flip() +
theme_bw() +
theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
ggarrange(p1, p2, nrow = 2, common.legend = TRUE, align = "v") %>%
annotate_figure(
top = text_grob("Distribución errores de validación cruzada", size = 15)
)El rmse (raíz del error medio cuadrático) promedio estimado mediante validación cruzada repetida es de 67990. Este valor será contrastado más adelante cuando se calcule el rmse del modelo con el conjunto de test.
Almacenar las predicciones de cada partición es útil para poder evaluar los residuos del modelo y diagnosticar así su comportamiento.
# Predicciones individuales de cada observación.
# Si summarize = TRUE se agregan todos los valores predichos a nivel de
# observación.
validation_fit %>% collect_predictions(summarize = TRUE) %>% head()p1 <- ggplot(
data = validation_fit %>% collect_predictions(summarize = TRUE),
aes(x = price, y = .pred)
) +
geom_point(alpha = 0.3) +
geom_abline(slope = 1, intercept = 0, color = "firebrick") +
labs(title = "Valor predicho vs valor real") +
theme_bw()
p2 <- ggplot(
data = validation_fit %>% collect_predictions(summarize = TRUE),
aes(x = .row, y = price - .pred)
) +
geom_point(alpha = 0.3) +
geom_hline(yintercept = 0, color = "firebrick") +
labs(title = "Residuos del modelo") +
theme_bw()
p3 <- ggplot(
data = validation_fit %>% collect_predictions(summarize = TRUE),
aes(x = price - .pred)
) +
geom_density() +
labs(title = "Distribución residuos del modelo") +
theme_bw()
p4 <- ggplot(
data = validation_fit %>% collect_predictions(summarize = TRUE),
aes(sample = price - .pred)
) +
geom_qq() +
geom_qq_line(color = "firebrick") +
labs(title = "Q-Q residuos del modelo") +
theme_bw()
ggarrange(plotlist = list(p1, p2, p3, p4)) %>%
annotate_figure(
top = text_grob("Distribución residuos", size = 15, face = "bold")
)Cuando se aplican métodos de* resampling debemos tener en cuenta dos cosas:
El coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto,
y la reproducibilidad en la creación de las particiones.
fit_resamples() permite paralelizar el proceso si se
registra un backed paralelo, para ello usamos la función
resterDoParallel() del paquete doParallel.
library("doParallel")Loading required package: foreach
Warning: package ‘foreach’ was built under R version 4.0.5
Attaching package: ‘foreach’
The following objects are masked from ‘package:purrr’:
accumulate, when
Loading required package: iterators
Loading required package: parallel
registerDoParallel(cores = detectCores() - 1)
set.seed(2020)
model_tree <- decision_tree(mode = "regression") %>%
set_engine(engine = "rpart")
validacion_fit <- fit_resamples(
object = model_tree,
preprocessor = transformer,
resamples = cv_folds,
metrics = metric_set(rmse, mae),
control = control_resamples(save_pred = TRUE)
)
stopImplicitCluster()Muchos modelos, entre ellos los árboles de regresión, contienen parámetros (datos de entrada al modelo) que no pueden aprenderse a partir de los datos de entrenamiento y, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. ¡Ojo! no se trata de tener un modelo y correrlo y ya con esto llamarnos científicos de datos.
Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades por medio de algoritmos de optimización.
Podemos seguir los siguientes pasos:
Escoger un conjunto de valores para el o los hiperparámetros.
Para cada valor (combinación de valores si hay más de un hiperparámetro), entrenar el modelo y estimar su error mediante un método de validación.
Ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.
El modelo decision_tree empleado hasta ahora tienen tres
hiperparámetros: cost_complexity, tree_depth y
min_n. Todos ellos controlan el tamaño del árbol
final. Por ejemplo, cuanto mayor el valor de
tree_depth, más ramificaciones tiene el árbol y más
flexible es el modelo. Esto le permite ajustarse mejor a las
observaciones de entrenamiento pero con el riesgo de
overfitting. Con las funciones tune() y
tune_grid() del paquete tune se pueden
explorar diferentes valores de hiperparámetros sin apenas tener que
cambiar la estructura del código.
Se vuelve a ajustar un modelo decision_tree con
diferentes valores de tree_depth y min_n, y
empleando validación cruzada repetida para identificar con cuáles de
ellos se obtienen mejores resultados.
# Definición del modelo y los hiperparámetros a optimizar
model_tree <- decision_tree(
mode = "regression",
tree_depth = tune(),
min_n = tune()
) %>%
set_engine(engine = "rpart")
# Validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = model_tree,
# El objeto recipe no tiene que estar entrenado
preprocessor = transformer,
# Las resamples se tienen que haber creado con los datos sin
# prerocesar
resamples = cv_folds,
metrics = metric_set(rmse, mae),
control = control_grid(save_pred = TRUE),
# Número de combinaciones generadas automáticamente
grid = 70
)
stopImplicitCluster()Los resultados de la búsqueda de hiperparámetros pueden verse
haciendo un unnest() de el tibble que se crea de la
ejecución grid_fit, o con las funciones auxiliares
collect_metrics(), collect_predictions,
show_best() y select_best().
grid_fit %>% unnest(.metrics) %>% head()grid_fit %>% collect_metrics(summarize = TRUE) %>% head()grid_fit %>% show_best(metric = "rmse", n = 5)Es conveniente visualizar la evolución del error en función de los hiperparámetros. Sin embargo, hay que tener en cuenta que estos no son independientes los unos de los otros, el comportamiento de un hiperparámetro puede cambiar mucho dependiendo del valor que tomen los otros.
grid_fit %>%
collect_metrics(summarize = TRUE) %>%
filter(.metric == "rmse") %>%
select(-c(.estimator, n)) %>%
pivot_longer(
cols = c(tree_depth, min_n),
values_to = "value",
names_to = "parameter"
) %>%
ggplot(aes(x = value, y = mean, color = parameter)) +
geom_point() +
geom_line() +
labs(title = "Evolución del error en función de los hiperparámetros") +
facet_wrap(facets = vars(parameter), nrow = 2, scales = "free") +
theme_bw() +
theme(legend.position = "none")grid_fit %>%
collect_metrics(summarize = TRUE) %>%
filter(.metric == "rmse") %>%
select(-c(.estimator, n)) %>%
ggplot(aes(x = tree_depth, y = min_n, color = mean, size = mean)) +
geom_point() +
scale_color_viridis_c() +
labs(title = "Evolución del error en función de los hiperparámetros") +
theme_bw()En este caso, parece que el error del modelo se reduce a medida que
el valor de min_n aumenta (aunque no es estable) y también
el de tree_depth. En este último la mejora parece
estabilizarse a partir de 4.
Aunque los autores de tidymodels han establecido por
defecto valores adecuados basados en su amplia experiencia, en muchos
casos es conveniente tener más control sobre la búsqueda, ya que la
selección se hace de forma automática. Para conseguir este control,
podemos utilizar expand.grid() (para especificar
exactamente qué valores se emplean), o con las funciónes
regular_grid(), grid_random(),
grid_max_entropy(), y grid_latin_hypercube()
del paquete dials.
Se repite la búsqueda de mejores hiperparámetros pero esta vez definiendo el espacio de búsqueda.
# Definición del modelo y de los hiperparámetros a optimizar
model_tree <- decision_tree(
mode = "regression",
tree_depth = tune(),
min_n = tune()
) %>%
set_engine(engine = "rpart")
# Validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Grid de Hiperparámetros
set.seed(1234)
hiperpar_grid <- grid_random(
# Rango de búsqueda para cada hiperparámetro
tree_depth(range = c(1, 10), trans = NULL),
min_n(range = c(2, 100), trans = NULL),
# Número combinaciones aleatorias probadas
size = 50
)
# Ejecución y optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = model_tree,
# El objeto recipe no tiene que estar entrenado
preprocessor = transformer,
# Las resamples se tienen que haber creado con los datos sin
# prerocesar
resamples = cv_folds,
metrics = metric_set(rmse, mae),
control = control_resamples(save_pred = TRUE),
# Hiperparámetros
grid = hiperpar_grid
)
stopImplicitCluster()Ahora vemos los resultados:
grid_fit %>% show_best(metric = "rmse", n = 5)grid_fit %>%
collect_metrics(summarize = TRUE) %>%
filter(.metric == "rmse") %>%
select(-c(.estimator, n)) %>%
pivot_longer(
cols = c(tree_depth, min_n),
values_to = "value",
names_to = "parameter"
) %>%
ggplot(aes(x = value, y = mean, color = parameter)) +
geom_point() +
geom_line() +
labs(title = "Evolución del error en función de los hiperparámetros") +
facet_wrap(facets = vars(parameter), nrow = 2, scales = "free") +
theme_bw() +
theme(legend.position = "none")Los mejores resultados, en base al rmse se han obtenido con los hiperparámetros:
# Selección de los mejores hiperparámetros encontrados
best_hip <- select_best(grid_fit, metric = "rmse")
best_hipYa haciendo la optimización de los hiperparámetros, entramos al modelo pero con los datos de entrenamiento con estos parámetros.
final_model_tree <- finalize_model(x = model_tree,
parameters = best_hip)
final_model_treeDecision Tree Model Specification (regression)
Main Arguments:
tree_depth = 10
min_n = 69
Computational engine: rpart
final_model_tree_fit <- final_model_tree %>%
fit(
formula = price ~ .,
data = data_train_prep
)Una vez que el modelo final ha sido ajustado, con la función
predict() se pueden predecir nuevos datos. Los argumentos
de la función son:
object: un modelo entrenado.
newdata: un dataframe con nuevas
observaciones.
type: tipo de predicción ("numeric",
"class", "prob", "conf_int",
"pred_int", "quantile", or
"raw").
predicciones <- final_model_tree_fit %>%
predict(
new_data = data_test_prep,
#new_data = bake(transformer_fit, datos_test),
type = "numeric"
)
predicciones %>% head()Ya hemos validado con modelos de validación (CV, bootstrapping…) y se
han conseguido estimaciones del error del modelo al predecir nuevas
observaciones. Pero la mejor forma de evaluar el modelo final es
prediciendo un conjunto test, es decir, utilizando un conjunto de datos
que no hemos usado en el entrenamiento y en la optimización. El paquete
yardstick
tiene funciones predefinidas para métricas en cuestión que nos pueden
ayudar.
predicciones <- predicciones %>%
bind_cols(data_test_prep %>% select(price))
rmse(predicciones, truth = price, estimate = .pred, na_rm = TRUE)En el apartado de optimización, se estimó, mediante validación cruzada repetida, que el rmse del modelo model_tree era de 67790, un valor próximo al obtenido con el conjunto de test baja a 61536.
workflows permiten combinar en un sólo objeto todos los
elementos que se encargan del procesamiento (recipes),
modelado (parsnip y tune) y post-procesado.
Para crear el workflows se van encadenando los elementos
con las funciones add_* o modificando los elementos de las
funciones update_*.
# Definición del modelo y los hiperparámetros a optimizar
model_tree <- decision_tree(
mode = "regression",
tree_depth = tune(),
min_n = tune()
) %>%
set_engine(engine = "rpart")
# Definición de procesado
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Workflow
workflow_modelado <- workflow() %>%
add_recipe(transformer) %>%
add_model(model_tree)
workflow_modelado══ Workflow ════════════════════════════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: decision_tree()
── Preprocessor ────────────────────────────────────────────────────────────────────────────────────────
5 Recipe Steps
• step_naomit()
• step_nzv()
• step_center()
• step_scale()
• step_dummy()
── Model ───────────────────────────────────────────────────────────────────────────────────────────────
Decision Tree Model Specification (regression)
Main Arguments:
tree_depth = tune()
min_n = tune()
Computational engine: rpart
El objeto workflow puede ajustarse directamente con la
función fit() o hacer tuning de hiperparámetros
con tune_grid().
# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
# Rango de búsqueda para cada hiperparámetro
tree_depth(range = c(1, 10), trans = NULL),
min_n(range = c(2, 100), trans = NULL),
# Número valores por hiperparámetro
levels = c(3, 3)
)
# Ejecución optimización hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = workflow_modelado,
resamples = cv_folds,
metrics = metric_set(rmse, mae),
control = control_resamples(save_pred = TRUE),
# Hiperparámetros
grid = hiperpar_grid
)
stopImplicitCluster()
# Entrenamiento final
best_hiper <- select_best(grid_fit, metric = "rmse")
final_model_fit <- finalize_workflow(
x = workflow_modelado,
parameters = best_hiper
) %>%
fit(
data = data_train
) %>%
pull_workflow_fit()Uno de los problemas que tenemos para encontrar los mejores hiperparámetros es el coste computacional. Una alternativa es la búsqueda de hiperparámetros con métodos de optimización bayesiana.
En términos generales, la optimización bayesiana de hiperparámetros consiste en crear un modelo probabilístico en el que la función objetivo es la métrica de validación del modelo (rmse, auc, precisión, etc). Con esta estrategia, se consigue que la búsqueda se vaya redirigiendo en cada iteración hacia las regiones de mayor interés.
El objetivo final es reducir el número de combinaciones de hiperparámetros con las que se evalúa el modelo, eligiendo únicamente los mejores candidatos. Esto significa que, la ventaja frente a las otras estrategias mencionadas, se maximiza cuando el espacio de búsqueda es muy amplio o la evaluación del modelo es muy lenta.
# Definición del modelo y los hiperparámetros a optimizar
model_tree <- decision_tree(
mode = "regression",
tree_depth = tune(),
min_n = tune()
) %>%
set_engine(engine = "rpart")
# Procesado
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# workflow
workflow_modelado <- workflow() %>%
add_recipe(transformer) %>%
add_model(model_tree)
# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
# Rango de búsqueda para cada hiperparámetro
tree_depth(range = c(1, 10), trans = NULL),
min_n(range = c(2, 100), trans = NULL),
# Número valores por hiperparámetro
levels = c(3, 3)
)
# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_bayes(
workflow_modelado,
resamples = cv_folds,
# Iniciación aleatoria con 20 candidatos
initial = 20,
# Numero de iteraciones de optimización
iter = 30,
# Métrica optimizada
metrics = metric_set(rmse),
control = control_bayes(no_improve = 20, verbose = FALSE)
)! No improvement for 20 iterations; returning current results.
stopImplicitCluster()Al igual que con el método anterior, podemos ver la evolución del error:
autoplot(grid_fit, type = "performance") +
labs(title = "Evolución del error") +
theme_bw()Para mostrar los mejores hiperparámetros encontrados
show_best(x = grid_fit, metric = "rmse")Nótese que incluso el rmse encontrado es más bajo que con el caso anterior. Por último realizamos el entrenamiento final:
best_hiper <- select_best(grid_fit, metric = "rmse")
final_model_fit <- finalize_workflow(
x = workflow_modelado,
parameters = best_hiper
) %>%
fit(
data = data_train
) %>%
pull_workflow_fit()Y evaluamos las predicciones
predicciones <- final_model_tree_fit %>%
predict(
new_data = data_test_prep,
#new_data = bake(transformer_fit, datos_test),
type = "numeric"
)
predicciones %>% head()Y por último encontramos el error de test:
predicciones <- predicciones %>%
bind_cols(data_test_prep %>% select(price))
rmse(predicciones, truth = price, estimate = .pred, na_rm = TRUE)Podemos llegar al mismo reslultado anterior de 61536 con un coste computacional menor.
En esta sección veremos modelos disponibles en la paquetería
parsnip y al final comparamos sus resultados, estos modelos
son:
GLM (Generalized Linear Models),
Random forest,
SVM (Suport Vector Machine), y
MARS (Multivariate adaptive regression splines).
Realizaremos el mismo procedimiento que hicimos con árboles de decisión, pero ahora cambiando al modelo de regresión lineal:
# Definición del modelo y los hiperparámetros (GLM)
model_glm <- linear_reg(
mode = "regression",
penalty = tune(),
mixture = tune()
) %>%
set_engine(engine = "glmnet", nlambda = 100)
# Procesado
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
# Validación y generación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Workflow
workflow_modelado <- workflow() %>%
add_recipe(transformer) %>%
add_model(model_glm)
# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
# Rango de búsqueda para cada hiperparámetro
penalty(range = c(0, 1), trans = NULL),
mixture(range = c(0, 1), trans = NULL),
# Número de combinaciones totales
levels = c(10, 10)
)
# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = workflow_modelado,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = TRUE),
# Hiperparámetros
grid = hiperpar_grid
)
stopImplicitCluster()Mostramos los mejores resultados de la optimización:
show_best(grid_fit, metric = "rmse", n = 10)Usamos los hiperparámetros óptimos para obtener el modelo final:
best_hiper <- select_best(grid_fit, metric = "rmse")
modelo_glm <- finalize_workflow(
x = workflow_modelado,
parameters = best_hiper
)
modelo_glm_fit <- modelo_glm %>%
fit(
data = data_train
)Calculamos las predicciones con el test
predicciones <- modelo_glm_fit %>%
predict(
new_data = data_test,
type = "numeric"
)Y por último calculamos el error
predicciones <- predicciones %>%
bind_cols(data_test_prep %>% select(price))
error_test_glm <- rmse(
data = predicciones,
truth = price,
estimate = .pred,
na_rm = TRUE
) %>%
mutate(
modelo = "GLM"
)
error_test_glmUn modelo Random Forest está formado por un conjunto de árboles de decisión individuales, cada uno ajustado empleando una muestra bootstrapping de los datos de entrenamiento.
En el entrenamiento de cada árbol, las observaciones se van distribuyendo por nodos generando la estructura del árbol hasta alcanzar un nodo terminal. Cuando se quiere predecir una nueva observación, esta recorre el árbol acorde al valor de sus predictores hasta alcanzar uno de los nodos terminales. La predicción del árbol es la media de la variable respuesta (la moda en problemas de clasificación) de todas las observaciones de entrenamiento que están en este mismo nodo terminal. La predicción de un modelo Random Forest es la media de las predicciones de todos los árboles que lo forman.
Encontremos el modelo utilizando random forest:
# Definición del modelo y los hiperparámetros a optimizar
modelo_rf <- rand_forest(
mode = "regression",
mtry = tune(),
trees = tune(),
min_n = tune()
) %>%
set_engine(engine = "ranger")
# Procesado
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Workflow
workflow_modelado <- workflow() %>%
add_recipe(transformer) %>%
add_model(modelo_rf)
# Grid de hiperparámetros
hiperpar_grid <- grid_max_entropy(
# Rango de búsqueda para cada hiperparámetro
mtry(range = c(1L, 10L), trans = NULL),
trees(range = c(500L, 3000L), trans = NULL),
min_n(range = c(2L, 100L), trans = NULL),
# Número de combinaciones totales
size = 100
)
# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = workflow_modelado,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = TRUE),
# Hiperparámetros
grid = hiperpar_grid
)
stopImplicitCluster()Mostramos los resultados de optimización:
show_best(grid_fit, metric = "rmse")Usamos los mejores hiperparámetros en el modelo:
best_hiper <- select_best(grid_fit, metric = "rmse")
modelo_rf <- finalize_workflow(
x = workflow_modelado,
parameters = best_hiper
)
modelo_rf_fit <- modelo_rf %>%
fit(
data = data_train
)Evaluamos las predicciones de test
predicciones <- modelo_rf_fit %>%
predict(
new_data = data_test,
type = "numeric"
)y encontramos el error
predicciones <- predicciones %>%
bind_cols(data_test_prep %>% select(price))
error_test_rf <- rmse(
data = predicciones,
truth = price,
estimate = .pred,
na_rm = TRUE
) %>%
mutate(
modelo = "RF"
)
error_test_rfEl método de clasificación-regresión Suport Vector Machine (SVM) es un método de no paramétrico, se fundamenta en el concepto de hiperplano y de máxima separación. Los modelos SVM encuentran el hiperplano óptimo capaz de predecir el valor/clasificación de las observaciones con el menor error posible.
Para que el modelo nos funcione correctamente debemos instalar el
paquete kernlab.
# Definición del modelo y los hiperparámetros a optimizar
modelo_svm <- svm_rbf(
mode = "regression",
cost = tune(),
rbf_sigma = tune(),
margin = tune()
) %>%
set_engine(engine = "kernlab")
# Procesado
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
# Estrategia de validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Workflow
workflow_modelado <- workflow() %>%
add_recipe(transformer) %>%
add_model(modelo_svm)
# Grid hiperparámetros
hiperpar_grid <- grid_random(
# Rango de búsqueda para cada hiperparámetro
cost(range = c(-10, -1), trans = log2_trans()),
rbf_sigma(range = c(-10, 0), trans = log10_trans()),
svm_margin(range = c(0, 0.2), trans = NULL),
# Número de combinaciones totales
size = 100
)
# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = workflow_modelado,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = TRUE),
# Hiperparámetros
grid = hiperpar_grid
)
stopImplicitCluster()Mostramos los resultados de la optimización
show_best(grid_fit, metric = "rmse", n = 10)Y utilizamos los hiperparámetros óptimos para el modelo
best_hiper <- select_best(grid_fit, metric = "rmse")
modelo_svm <- finalize_workflow(
x = workflow_modelado,
parameters = best_hiper
)
modelo_svm_fit <- modelo_svm %>%
fit(
data = data_train
)Encontramos las predicciones en el test:
predicciones <- modelo_svm_fit %>%
predict(
new_data = data_test,
type = "numeric"
)Y finalmente encontramos el error
predicciones <- predicciones %>%
bind_cols(data_test_prep %>% select(price))
error_test_svm <- rmse(
data = predicciones,
truth = price,
estimate = .pred,
na_rm = TRUE
) %>%
mutate(
modelo = "SVM"
)
error_test_svmMultivariate adaptive regression splines (MARS) es una extensión no paramétrica de los modelos de regresión lineal que automatiza la incorporación de relaciones no lineales entre los predictores y la variable respuesta, así como interacciones entre los predictores.
Para que este funcione bien debemos instalar el paquete
earth.
# Definición del modelo y los hiperparámetros a optimizar
modelo_mars <- mars(
mode = "regression",
num_terms = tune(),
prod_degree = 2,
prune_method = "none"
) %>%
set_engine(engine = "earth")
# Procesado
transformer <- recipe(
formula = price ~ .,
data = data_train
) %>%
step_naomit(all_predictors()) %>%
step_nzv(all_predictors()) %>%
step_center(all_numeric(), -all_outcomes()) %>%
step_scale(all_numeric(), -all_outcomes()) %>%
step_dummy(all_nominal(), -all_outcomes())
# Validación y creación de particiones
set.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
strata = price
)
# Workflow
workflow_modelado <- workflow() %>%
add_recipe(transformer) %>%
add_model(modelo_mars)
# Grid de hiperparámetros
hiperpar_grid <- grid_regular(
# Rango de búsqueda para cada hiperparámetro
num_terms(range = c(1, 20), trans = NULL),
levels = 20
)
# Optimización de hiperparámetros
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
object = workflow_modelado,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = TRUE),
# Hiperparámetros
grid = hiperpar_grid
)
stopImplicitCluster()Mostramos los resultados de la optimización
show_best(grid_fit, metric = "rmse", n = 10)Utilizamos en el modelo los hiperparámetros óptimos:
best_hiper <- select_best(grid_fit, metric = "rmse")
modelo_mars <- finalize_workflow(
x = workflow_modelado,
parameters = best_hiper
)
mmodelo_mars_fit <- modelo_svm %>%
fit(
data = data_train
)Encontramos las predicciones test:
predicciones <- mmodelo_mars_fit %>%
predict(
new_data = data_test,
type = "numeric"
)Encontramos el error
predicciones <- predicciones %>%
bind_cols(data_test_prep %>% select(price))
error_test_mars <- rmse(
data = predicciones,
truth = price,
estimate = .pred,
na_rm = TRUE
) %>%
mutate(
modelo = "MARS"
)
error_test_marsset.seed(1234)
cv_folds <- vfold_cv(
data = data_train,
v = 5,
repeats = 5,
strata = price
)
registerDoParallel(cores = parallel::detectCores() - 1)
validacion_glm <- fit_resamples(
object = modelo_glm,
# preprocessor = transformer,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = FALSE)
) %>%
collect_metrics(summarize = FALSE) %>%
mutate(modelo = "GLM")
validacion_svm <- fit_resamples(
object = modelo_svm,
#preprocessor = transformer,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = FALSE)
) %>%
collect_metrics(summarize = FALSE) %>%
mutate(modelo = "SVM")
validacion_rf <- fit_resamples(
object = modelo_rf,
# preprocessor = transformer,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = FALSE)
) %>%
collect_metrics(summarize = FALSE) %>%
mutate(modelo = "RF")
validacion_mars <- fit_resamples(
object = modelo_mars,
# preprocessor = transformer,
resamples = cv_folds,
metrics = metric_set(rmse),
control = control_resamples(save_pred = FALSE)
) %>%
collect_metrics(summarize = FALSE) %>%
mutate(modelo = "MARS")
stopImplicitCluster()bind_rows(
validacion_glm,
validacion_rf,
validacion_svm,
validacion_mars
) %>%
ggplot(aes(x = modelo, y = .estimate, color = modelo)) +
geom_violin() +
geom_boxplot(outlier.shape = NA, width = 0.2) +
geom_point(alpHa = 0.1) +
labs(title = "Error de validación cruzada", y = "RMSE") +
theme_bw() +
theme(legend.position = "none")Warning: Ignoring unknown parameters: alpHa
errores_test <- bind_rows(
error_test_glm,
error_test_rf,
error_test_svm,
error_test_mars
)
errores_test %>% select(modelo, .metric, .estimate)ggplot(data = errores_test) +
geom_col(aes(x = modelo, y = .estimate, fill= modelo), color = "gray") +
coord_flip() +
labs(title = "Error de test") +
theme_bw()Kuhn et al., (2020). Tidymodels: a collection of packages for modeling and machine learning using tidyverse principles. https://www.tidymodels.org
G. James, D. Witten, T. Hastie, R. Tibshirani. An Introduction to Statistical Learning. MIT Press, 2010.
Linear Models with R By Julian J. Faraway
Machine Learning con R y tidymodels, Joaquín Amat Rodrigo.